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Abstract 

A 3D simulation of a non-relativistic, magnetically driven jet propagating in a stratified atmosphere is presented, covering about three 
decades in distance and two decades in sideways expansion. The simulation captures the jet acceleration through the critical surfaces 
and the development of (kink-)instabilities driven by the free energy in the toroidal magnetic field component. The instabilities destroy 
the ordered helical structure of the magnetic field, dissipating the toroidal field energy on a length scale of about 2-15 times the Alfven 
distance. We compare the results with a 2.5D (axisymmetric) simulation, which does not become unstable. The acceleration of the 
flow is found to be quite similar in both cases, but the mechanisms of acceleration differ. In the 2.5D case approximately 20% of the 
Poynting flux remains in the flow, in the 3D case this fraction is largely dissipated internally. Half of the dissipated energy is available 
for light emission; the resulting radiation would produce structures resembling those seen in protostellar jets. 
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1. Introduction and rationale of the calculations 

A magnetized outflow produced by a rotating magnetic object 
has become the default interpretation for objects ranging from 
protostellar jets to gamma-ray bursts. The outflow in this model 
contains a tightly wound helical magnetic field. Such a nearly 
toroidal field represents a source of free energy that makes the 
flow inherently prone to non-axisymmetric magnetic instabili- 
ties. In this study we investigate the longer-term development of 
such instabilities and their consequences for jet phenomenology. 

Instabilities are not necessarily fatal for the jet. Kink insta- 
bilities of helical magnetic field configurations typically saturate 
at a finite amplitude. Such instabilities at moderate amplitudes 
have been invoked to explain phenomena like the wiggly appear- 
ance of Herbig-Haro objects (Todo et al. 1993) or the orientation 
of VLBI jets in AGN (Konigl & Choudhuri 1985). 

The development of kink instability in a jet is not the same 
as in a laboratory configuration. Due to the sideways expansion 
of the flow, the ratio of poloidal (stabilizing) to toroidal (destabi- 
lizing) field components decreases with distance along the flow. 
Conditions favorable for instability are thus continually recre- 
ated in such a flow. In Moll et al. (2008), hereafter Paper 1, we 
presented 3D MHD simulations showing the onset of instabili- 
ties in an expanding jet created by twisting a purely radial mag- 
netic field. The degree of instability was found to depend on the 
kind of rotation that generates the twist. The highest degree of 
instability was attained with a constant angular velocity (rigid 
rotation); the jet produced in this case was subject to helical 
deformations with large amplitudes, causing sideways displace- 
ments of several degrees. However, the magnetic structure of the 
jet was not disrupted within the computational volume, and the 
instability did not lead to a significant decrease of the Poynting 
flux. These results indicate the need to follow the instabilities to 
larger distances from the source. We present here the results of 
simulations extending to a distance of 1000 times the diameter 
of the jet source. 
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In Paper 1 we showed how the degree of instability depends 
on the way in which the jet is collimated by its environment. 
If collimation conditions are such that the opening angle of the 
jet increases with distance, the Alfven travel time across the jet 
increases more rapidly than the expansion time scale. As a re- 
sult, instabilities soon "freeze out", and decay of the toroidal 
field becomes ineffective. In better collimating jets, such that the 
opening angle narrows with distance, instability is always effec- 
tive. The calculations presented here are for such a case. It is 
probably the most relevant for both AGN and protostellar jets. 
Observations of the jets in M87 (Junor et al. 1999) and HH30 
(Mundt et al. 1990), for example, show a rapidly decreasing 
opening angle in the inner regions of the jet. 

The dissipation of the magnetic energy in the toroidal field 
heats the plasma. Radiation produced by this plasma would be an 
alternative or complement to the standard mechanism invoked, 
which relies on dissipation in internal shocks. Observations 
such as those of M87 and HH30 indicate that dissipation starts 
fairly close to the central object, compared with most observ- 
able length scales, but also that the innermost regions, perhaps 
comparable to the Alfven distance r A , are quiet. Inferences from 
VLBI observations indicate that AGN jets are not magnetically 
dominated on most observable scales (cf. Sikora et al. 2005, and 
references therein); decay of the magnetic field relatively close 
to the source of the jet would fit this observation (see also dis- 
cussion in Giannios & Spruit 2006). 

If dissipation takes place close enough to the source, it is pos- 
sible that it can be studied realistically by 3D simulations with 
current computational resources. In Paper 1, we already found 
indications that decay of the toroidal field may become impor- 
tant at distances as close as 10-30 times r A . 

The release of magnetic energy by the instability may 
also be important for accelerating the flow (Drenkhahn 2002). 
Dissipation of toroidal field causes the magnetic pressure gradi- 
ent along the jet to steepen; this adds an accelerating force that 
is absent when the toroidal field is conserved (Giannios & Spruit 
2006; Spruit 2008). 
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To facilitate extraction of physics from the numerical results, 
the 3D results in the following are compared systematically with 
a 2.5D (axisymmetric) simulation corresponding to the same ini- 
tial and boundary conditions. 



2. Methods 

2.1. Numerical MHD solver, grid and coordinates 

The following is a brief summary; for details see Paper 1 where 
the same numerical approach but different initial conditions 
were used. Details on the MHD code can also be found in 
Obergaulinger (2008). 

We numerically solve the ideal adiabatic MHD equations in 
a static external gravitational potential (D oc r _1 on a spherical 
grid (r, 9, <p). In the 2.5D, axisymmetric simulation, the jet prop- 
agates along the coordinate axis = 0. In the 3D simulation, 
the jet's axis is taken in the direction 6 = (f> = n/2. The jet thus 
propagates in equatorial direction of the coordinate system. This 
avoids the coordinate singularity along the polar axis (9 = 0), 
which is numerically problematic for trans-axial flows such as 
are caused by instabilities. The computational volume covers a 
range A9 = A(p that comprises about twice the expected open- 
ing angle of the jet. The spacing of the computational grid is 
uniform in the angular directions and logarithmic in the radial 
direction. In this way, the varying numerical resolution approx- 
imately matches the increase of the natural length scales in the 
expanding jet. 

For presentation and discussion, we transform the results to a 
different coordinate system. This is again a spherical coordinate 
system, but with the polar axis aligned with the jet. The polar 
and azimuthal angles in this system are denoted by § and <p, 
respectively: 



sin 2 -& = 1 - sin 2 9 sin 2 (p, 
tan <p = tan 9 cos (p 



(1) 
(2) 



in the 3D simulation and -& = 9, <p = <p in the 2.5D simulation. 
R .-r sin § is used to denote the distance to the axis (cylindrical 
radius). 

Dissipation of magnetic energy in the flow causes heating. 
In nature this would lead to losses by radiation; in the compu- 
tations it can cause numerical problems in regions where the 
magnetic energy dominates. Instead of a more realistic model 
for such losses, a temperature-control term is added in the en- 
ergy equation, such that the temperature relaxes to that of the 
initial state on an appropriate time scale. In the simulations pre- 
sented here this time scale is chosen such that the temperature 
stays within about a factor 100 around the initial value. 

2.2. Initial and boundary conditions 

The initial state consists of a current-free magnetic field embed- 
ded in a stratified atmosphere in the gravitational field of a point 
mass. The field configuration of this initial state is of a "collimat- 
ing" type, the distance between neighboring field lines increases 
less rapidly than the distance from the source r. In Paper 1, we 
showed that instabilities develop more strongly in such collimat- 
ing configurations than in a purely radial initial field. 



2.2.1. Initial field configuration 

The initial field is axisymmetric around the jet axis and hence 
can be written as 

fi = ^x^ = Vx(^)=VxA, (3) 



R 



where if/ is the stream function of the field. We construct the 
initial condition as a linear combination of the stream function 
of a monopole field, given by 



lAmono K 1 - COS 

and a field with the stream function 



lApara <* 




(4) 



(5) 



(z > 0) of which the field lines have a parabolic shape (Cao & 
Spruit 1994). Here, z = r cos & denotes the height along the jet's 
central axis (which is actually the y-axis in the 3D simulation 
and the z-axis in the 2.5D simulation). The weighting of the two 
components is parametrized with the quantity 



Br, 



(6) 



r=rb,i?=0 



the relative strength of the constituent fields at the lower bound- 
ary r = r\,. £ — > oo corresponds to a pure monopole field and 
£ = corresponds to a pure parabolic field. The radius r equ at 



which B mono = B para on the central axis follows from 



r equ(^0 + r n ) 



rl(Ro + r equ ) 



(7) 



The lower the value of £ or r e q U , the more collimating is the shape 
of the initial field. 



2.2.2. Stratification 

We impose the static gravitational field of a point mass, located 
at the origin of the coordinate system. The stratification is ini- 
tially in hydrostatic equilibrium in this potential. The gas pres- 
sure is chosen such that the plasma-/? := p/p m ag in the initial 
state is approximately constant at small radii (<«cr equ ), where the 
monopole field dominates: p oc r~ 4 , p oc r~ 3 , T oc r~ l , c s oc r -1 ^ 2 
and va k r~ 1/2 in the initial state. At large radii, the value of the 
initial ft is reduced as the magnetic field strength decreases less 
rapidly than r~ 2 . 



2.2.3. Boundaries 

Boundary conditions are maintained through the use of "ghost 
cells" outside the computational domain. At the sides (9 in the 
2.5D simulation, 9 and <f> in the 3D simulation) and top (upper r), 
we use open boundary conditions that allow for an almost force- 
free outflow or inflow of material, including magnetic fields. 

The bottom boundary is located at a finite height r\, above 
the origin of the gravitational potential. This distance is 1/200 
of the size of the computational volume. The jet is generated 
there by a "rotating disk" 1 of radius R\, around the axis. It is 



1 Strictly speaking, because of the spherical grid, what we call "disk" 
here is really a spherical cap with small curvature. 
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Table 1. Normalization units 
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implemented by maintaining a velocity field v = v»$n in the 
ghost zones: v v = v™™R/Rb for R < R\, and elsewhere. The 
disk thus rotates rigidly like in the cases R2 and R3 in Paper 1 . 

All quantities at the bottom boundary except for B are fixed 
at their initial values in the ghost cells. B is extrapolated from 
the interior of the domain. 



2.3. Parameters and units 

The above suggests a specification of the problem in terms of 7 
parameters: the field configuration parameter £ a scale for the 
field strength, scales for the pressure and density, the bottom 
boundary location rt,, the radius R\, of the rotating disk, and its 
rotation rate. Because of the symmetries of the problem, it is ac- 
tually determined by only 4 parameters; the remaining 3 depen- 
dences are equivalent to scaling factors. As the 4 independent 
parameters we choose the following dimensionless quantities: 
the plasma-/? of the initial state, the angle 0b = arctan(/?b/rb) 
which controls the opening angle of the flow, the field configu- 
ration parameter and finally the Alfvenic Mach number Ma of 
the rotating velocity field at the edge (R = R\,) of the launching 
disk, which controls the power of the jet. They have the values 
P = 1/9, £ = 30, 0b = 5.7°, M A = 0.1. [The simulations re- 
ported in Paper 1 corresponded to /? = 1/9, f — > oo, b = 5.7°, 
M A = 0.1.] 

The units used for reporting the results below are the length 
scale Iq = 2Rb, and the pressure and density on the axis at the 
bottom boundary, po = pb, po = pb- Together with the 4 model 
parameters p\ §b an d M A , the units for other quantities follow 
from these as shown in Table 1. 

The simulations cover a distance of 2000 times the initial jet 
radius Rb in the spherical range 5 < r < 1005. The resolution 
used in the 3D simulation is 768 x 128 x 128; the correspond- 
ing domain size is 1000 x 33.8° x 33.8°. The resolution used in 
the 2.5D simulation is 768 x 96; the corresponding domain size 
1000 x 16.8°. In both simulations, the radial width Ar of the grid 
cells increases from 0.03 at the lower boundary to 6.92 at the 
upper boundary in logarithmic steps. 



3. Results 

The 3D simulation was done using MPI parallelization on 128 
CPUs, the 2.5D simulation was done with OpenMP paralleliza- 
tion on 32 CPUs. 

The jet crossed the upper boundary of the computational vol- 
ume at the physical time t m 523 on the 13th wall clock day of 
the 3D simulation. We stopped it after 26 days, at which time 
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Figure 2. Volume rendered image of the 3D jet that shows what it might 
look like in observations. For the volume rendering, a simple model was 
used in which emissivity and opacity depend on temperature and mag- 
netic field strength. The jet exhibits a wiggly structure with bright knots, 
produced by the instabilities and the dissipation of magnetic energy. 



t = 1055 had been reached. The 2.5D simulation ran for 24 
hours, reaching t = 1732. 

The 3D jet is subject to non-axisymmetric instabilities, ev- 
idently of the kink (m = 1) kind. They have a disruptive effect 
on the magnetic field structure and cause the toroidal field to 
decay, see Fig. 1. If such a jet was observed, it would proba- 
bly look similar to Fig. 2, with bright knots and wiggles being 
prominent features. The knots move at a substantial fraction of 
the flow speed, sometimes merging or fading before leaving the 
computational domain. The 2.5D jet does not exhibit any form 
of instability. 

3.1. Acceleration, collimation and mass flow 

Most of the acceleration takes place below r = 100, where the 
jet reaches >80% of its terminal speed, see Fig. 3. The loca- 
tion of the sonic (v r — c s ~ c s v A j(cl + v A )'^ 2 , the slow mag- 



netosonic cusp speed), Alfven (v r = v Ar ) and fast magnetosonic 
(vj. = c 2 s + v 2 ) radii depends on the direction 0, see inset in Fig. 3. 
This is because the acceleration is more effective near the bound- 
ary of the jet, which is a consequence of the rigid rotation profile 
used in its generation, and because the poloidal magnetic field is 
being redistributed such that the local Alfven velocity becomes 
relatively high near the axis. As the toroidal field energy, which 
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Figure 1. a) Selected magnetic field lines in the 3D simulation on successively increasing length scales. The color coding gives the magnetic field 
strength relative to its initial value. The field lines shown are the ones that are anchored in the rotating disk at the lower boundary. The jet starts 
out with a helical magnetic field {left image) whose toroidal component becomes increasingly stronger (middle image) until instabilities disrupt 
the ordered structure, the toroidal field decays and the field becomes predominantly poloidal (right image), b) Radial component of the current 
density (V x B) in the 3D simulation. The backward current, shown in blue here, is concentrated on the central axis until being first displaced and 
then disrupted by instabilities. 
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Figure 3. Maximum velocity in jet direction as a function of distance. 
The inset shows the location of the critical surfaces. The flow passes first 
the sonic, then the Alfven and finally the fast magnetosonic surface. 



determines the jet acceleration and the development of instabili- 
ties, is concentrated towards the jet boundary (as opposed to the 
axis), the Alfven radius there is probably the most important for 
the subsequent considerations. 

The central velocity tends to be higher in the 3D simulation, 
presumably due to a more effective transfer of momentum from 
the boundary to the center. This may, together with the entrain- 
ment of ambient material discussed below, also explain why the 
peak velocities are somewhat lower in the 3D simulation. 
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Figure 4. Jet boundary, defined by the angle that encloses 95% of the 
total energy flow (see Eq. 8), as a function of distance. The unstable 3D 
jet is less collimated than the 2.5D jet. 



The instabilities have a noticeable effect on the collimation 
behavior, see Fig. 4. In the 2.5D simulation, the opening angle of 
the jet decreases by about 2° in the first half of the computational 
volume, which is about 1° less than what is marked by the shape 
of the initial magnetic field. In the second half, well beyond the 
Alfven surface, the opening angle settles to a constant value. The 
jet in the 3D simulation is less collimated 2 , the location of the 
boundary fluctuates with time in the unstable region. Averaged 
over time, we find it to be nearly conical. 

The mass flow rate M(t, r) := f pv r dA is somewhat 
higher in the 3D simulation and subject to strong fluctuations 
above r ~ 200: the mean value over all radii at t » 1060, mea- 
sured in the unit for M listed in Table 1, is 0.082 ± 34% in the 
2.5D and 0.082 + 34% in the 3D simulation. The average over 



2 Note, however, that with our definition of the jet boundary, a rigid 
displacement away from the center also increases the opening angle. 
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Figure 5. Energy flow rates through the r = const, surface. Magnetic 
enthalpy (integral over Poynting flux, blue line) is reduced in the 3D 
simulation as B 1 decays as a consequence of instabilities. 



several time steps shows a slight increase of M with r in the 
3D case, whereas no trend can be deduced in the 2.5D case. This 
may be an indicator for an enhanced entrainment of ambient ma- 
terial caused by the instabilities. It may also in part explain the 
lower peak velocities in the 3D simulation. 



3.2. Energy 

A comparison of the different kinds of energy flow rates gives 
information about energy transformations taking place in the jet. 
Integrating the radial energy flux 



+ p<Dv,. + S r (8) 

grav. potential magnetic enthalpy 
thermal enthalpy 




kinetic 



over the r — const, surface, we obtain the rate of energy fitotM 
that crosses the surface. The individual components, which we 
denote with £kj n , £thrm, Sgrav and £ mag , are plotted in Fig. 5. As 
the values are strongly fluctuating in the unstable region, we av- 
eraged over time in the 3D case. £ mag is reduced by 80% in the 
2.5D simulation and by 94% in the 3D simulation along the sim- 
ulated distance, the gap between the curves widens considerably 
where the 3D jet exhibits instabilities. We presume this to be 
caused by (numerical) magnetic dissipation, which turns mag- 
netic energy into heat, causing an increase in £ t hr m . The environ- 
mental thermal energy is not affected much. As our calculations 
include an energy loss term (see Sect. 2. 1), there is no one-to-one 
correspondence between the dissipated energy and the increase 

in fithrm- 

The final contribution of <5 mag to <5 tot is 24% in the 2.5D sim- 
ulation and 9% in the 3D simulation. S^n is up to about 10% 
smaller in the 3D simulation and does not show a dissipation- 
induced increase. arises from azimuthal and radial motion, 
with the relative share being highly similar in the 2.5D and 3D 
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Figure 6. Red: magnetic flux S(r) within 5.7° from the jet axis, nor- 
malized by its initial value. Green: flux of B* > (B = fi* + B~), see 
Sect. 3.3. 



case: the azimuthal contribution drops continuously to values 
<5% at r > 200 and the contribution from v§ is insignificant 
throughout. 

£ mag can be decomposed further. The radial component of 
the Poynting vector has 4 terms, 

S r = ±- (Blv r + B\v r - B § B r v § - B^B rVip ) . 



(9) 



Integrating these terms over the r = const, surface gives the com- 
ponents [in order of their appearance in Eq. (9)] . . . £ mao 
which, added together, give £ mag . These components of the mag- 
netic enthalpy flow rate are plotted in the lower panel of Fig. 5. 
At very small radii, the most important term is £ mag , the work 
done by the azimuthal flow against the azimuthal component of 
magnetic stress. Its contribution to £ mag decreases rapidly, how- 
ever, falling below 50% at r « 30 (i.e. well below the Alfven 
surface, near the sonic one). It is then taken over by £ mac „ which 
describes the flow of magnetic enthalpy stored in the azimuthal 
field. fi mag has minor significance in the 3D simulation, which 
may be due to the perturbed toroidal field having also a B§ com- 
ponent (a rigid displacement of a pure azimuthal field in the 
r = const, plane introduces a non-azimuthal component). £ mao 
is insignificant in both cases. The strong decrease of £ mag in the 
3D simulation is caused by the decrease of £ ma „. 

We find a net outflow of magnetic enthalpy through the lat- 
eral boundaries of the computational volume at the height of the 
jet front, with peak rates of about 0.15 in the 2.5D case and 0.05 
in the 3D case. The outflow is transient in the 2.5D case, vanish- 
ing quickly after the jet front leaves the computational volume. 
In the 3D case, however, it turns into an inflow of the order -0.05 
which persists until the end of the simulation. The energy in the 
radial magnetic field increases correspondingly, mainly outside 
the jet at# > 5.7°. 



3.3. Magnetic field: poloidal vs. toroidal 

The red line in Fig. 6 shows the magnetic flux S contained within 
an angle i? < 5.7° from the axis, 



E(r,t) 



X 

^ .'J 



const. 

»<5.T 



B r (f)dA, 



(10) 



divided by its initial value. For comparison, the green curve 
shows H + , the flux of only those field lines that have the same 
direction as the initial field (the green and the red curves co- 
incide in the 2.5D case). Although the 2.5D jet fills only part of 
the 5.7° cone, it causes a significant reduction of S . This reflects 
the lateral expansion of the jet due to the pressure exerted by B^. 
There is less reduction in the 3D case, with E being near the 
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Figure 9. Poloidal magnetic field lines and radial velocity 
field, plotted as a function of direction & (horizontal axis) 
and distance r (vertical axis); a conical magnetic surface 
would appear as a vertical line. The left-hand plot shows 
a snapshot of the 2.5D case, the right-hand plot shows the 
3D case with the magnetic field and velocity being aver- 
aged in azimuthal direction as well as time. Both plots 
show the same field lines in the sense that they start from 
the same direction at the lower boundary. The angular 
separation between field lines in the outer, high-?? part of 
the jet (here: between the 4th and 5th line from the left) 
increases noticeably in the 2.5D case, leading to acceler- 
ation by the magnetic "diverging nozzle" effect there, see 
discussion. 



dotted: 2^D @ 806<t<1059 solid: 3D @ 803<t<1055 



dotted: 2aD @ 806<t<1059 solid: 3D @ 803<t<1055 



5 r 

4 j- 

3 r 

2 r 



\ 



D / B 2 (t = 0,iS = 0) 

(Bp) / B 2 (t = 0,i? = 0) 





1000 



Figure 8. Magnetic field components, integrated over the jet cross sec- 
200 400 600 800 1 000 tion, as a function of distance, on a logarithmic scale. In a spherically 

r expanding ballistic flow these quantities would be constant, see text. 



Figure 7. Mean magnetic energies in the poloidal and toroidal fields 
within the jet, normalized by the initial on-axis magnetic energy, and 
their ratio. Without instability (2.5D), the field becomes predominantly 
toroidal at large distances. With instability (3D), it becomes predomi- 
nantly poloidal. 

original value at very large distances. The difference becomes 
evident where the 3D jet is unstable, showing that it is due to the 
dissipation of the toroidal field. There is only a small difference 
between the red and green curve in the 3D case, meaning that 
negative values of B r contribute little to the net flux. This shows 
that the toroidal field component dissipates without producing 
much "tangling" of the field lines. 

Fig. 7 compares the mean energy in the poloidal and toroidal 
magnetic fields over the width of the jet as a function of dis- 
tance, the width of the jet being defined by a suitable veloc- 
ity threshold. The 3D and 2.5D jets start out similar, with the 
mean poloidal and toroidal fields becoming comparable near the 
Alfven surface. Beyond that distance, the mean toroidal field en- 
ergy increases more strongly in the 2.5D case, roughly propor- 
tional to r~ 2 . In the 3D case, the instability-induced destruction 
of the toroidal field causes the slope to steepen substantially at 
r > 200, approximately as (B 2 ) ~ r~ 3 . The mean poloidal field 
energy, on the other hand, stays near the initial value in the 3D 



case while decreasing somewhat in the 2.5D case. Taking the 
ratio between the energies, we find the magnetic field to be pre- 
dominantly toroidal in the 2.5D case and predominantly poloidal 
in the 3D case at large distances. 

We also experimented with means of the form 
JXv r dA/ J v r dA, where the integral is performed over 
the whole r = const, surface, and obtained similar results. 
Taking the mean within the static # < 5.7° cone instead of using 
a velocity threshold yields a smaller value for (B£)/(B?) in the 
2.5D case, viz. 2 + 0.2 for r > 200. This is because an angle of 
5.7° includes more of the environment of the jet (cf. Fig. 4). 

In a ballistically expanding jet (constant velocity and open- 
ing angle) the magnetic field components vary with distance r as 
B r ~ r~ 2 , Be ~ B,p ~ r . Integrals over the width of the jet of 
B r , B^ and B 2 g are then constants. These integrals are shown in 
Fig. 8. The ballistic approximation works well in the 2.5D case, 
for B r as well as for B v , if the acceleration region is excluded. In 
the 3D case however, the integral of B r increases by about an or- 
der of magnitude along the jet, while the integral of B 2 decreases 
over the entire range. 

Fig. 9 shows some poloidal magnetic field lines and the ra- 
dial velocity in the 2.5D and 3D jet. The latter was "axisym- 
metrized" by averaging over the azimuthal coordinate <p. In the 
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Figure 10. Rate of work done in the volume delimited by the lower 
boundary and r by the components of the Lorentz force in the direction 
of the flow. The power delivered by the total Lorentz force is about the 
same in the 3D and 2.5D cases. 



outer, high-)?, part of the 2.5D jet, where the toroidal field is 
especially strong and jet acceleration most efficient, there is an 
increase of the angular separation between the field lines which 
is absent in the 3D jet. In other words, the poloidal magnetic flux 
decreases locally faster with distance in the 2.5D case. 

3.4. Forces and powers 

To identify the accelerating forces that generate the kinetic en- 
ergy flow discussed in Sect. 3.2, we compute P — J F v dV, the 
instantaneous power (rate of work) delivered by a specific force 
F in the direction of the flow in the integrated volume. The com- 
bination of gas pressure and gravitational forces, -Vp - pVcD, 
accounts for about one third of the power delivered by the sum 
of all forces in the whole volume, the corresponding accelera- 
tion takes place mainly below r ~ 30 (sonic surface). The rest 
is accounted for by the Lorentz force, which we decompose as 
follows: 



(V x B p ) x B p + (V x B v ) x 



+(V x B p ) xB^ + iVx B v ) x B p 



(11) 



with B^ = B^e^ and B p 



B - B^ The corresponding com- 
ponents of P are plotted in Fig. 10 as a function of the upper 
integral limit. 

The last force in Eq. (11) has only an azimuthal component. 
It is important mainly below the Alfven radius (r ~ 60), exert- 
ing a torque in the same direction in which rotation is applied 
at the lower boundary. The next-to-last force vanishes in the ax- 
isymmetric case. In the general case, it has only non-azimuthal 
components. Unlike the last force, it works against the flow; 
the two forces largely cancel each other in the 3D simulation. 
The Lorentz force associated with [second term in Eq. (11)] 
performs the same work in the 3D and 2.5D cases, despite the 
steepening of the toroidal magnetic pressure profile caused by 
the dissipation. The total power is also similar in the two cases, 
in agreement with the similarity of the kinetic energy flows in 
Fig. 5. 



The situation is a bit different if only the radial power 
J F, v r dV is taken into account. The Lorentz force associated 
with By does approximately 20% more work in the whole vol- 
ume in the 3D case, the additional power is delivered above 
r > 200. The power of the net force, however, is approximately 
the same in both cases. 



4. Summary and discussion 

We have simulated jets generated by twisting a parabolically 
shaped large-scale magnetic field in both 3D and axisymmetric 
2.5D. The shape of the initial field reflects itself in a fair amount 
of jet collimation, with an opening angle that decreases with dis- 
tance, thus facilitating the growth of instabilities. The simula- 
tions cover the acceleration phase, where the jet passes through 
the critical surfaces (sonic, Alfven and fast magnetosonic) for 
stationary MHD flows, as well as a substantial distance beyond 
these. We thus observe the onset of kink instabilities above the 
Alfven surface in the 3D simulation. The instabilities disrupt the 
magnetic field structure. They cause magnetic dissipation, sig- 
nificantly reducing the toroidal field strength and with it the flow 
rate of magnetic enthalpy (surface integral of the Poynting flux) 
on a length scale of about 2-15 times the minimal Alfven dis- 
tance in the jet. 

A direct comparison of the 2.5D and 3D simulations reveals 
no significant difference in the way the kinetic energy of the 
jet (integrated over its cross section) increases with distance. 
However, the distribution of the kinetic energy across the jet in- 
dicates differences in the acceleration process (cf. Fig. 9). In the 
axisymmetric flow, the acceleration is restricted to magnetic sur- 
faces that diverge from each other more rapidly than in a flow 
with fixed opening angles. This creates a "magnetic nozzle" ef- 
fect (Begelman & Li 1994) restricted to a limited range of angles 
within the flow (cf. discussion in Spruit 2008). A similar case of 
non-uniform expansion has been found in relativistic flows by 
Tchekhovskoy et al. (2008). In the 3D case, the (time averaged) 
acceleration is more uniform across the jet. It thus seems that 
while the steepening of the magnetic pressure gradient caused 
by the dissipation of the toroidal field leads to additional accel- 
eration in the 3D case, this is made up for by a more favorably 
distributed poloidal magnetic flux in the 2.5D case which also 
enhances acceleration. In other words, there are different mech- 
anism at work which yield the same result. Whether it is coin- 
cidence that the two effects have nearly the same strength is not 
clear. 

The energy released by the dissipation of the toroidal mag- 
netic field may be radiated as light. The yield depends on the de- 
tails of the dissipation (reconnection) and radiation mechanisms 
involved. Since the magnetic energy density (B^l%n) accounts 
for half of the Poynting flux (ByV r /4n), the present simulations 
suggest that the available luminosity may be as much as 10% of 
the initial magnetic enthalpy flow rate (Fig. 5), provided that half 
of it is converted into kinetic energy and the rest into light. As 
dissipation is not a smooth process, the emission will be stronger 
in some regions. These would likely turn up as bright knots in 
observations. Together with the wiggles caused by the kink in- 
stabilities, the knots produce a structured appearance that is sim- 
ilar to what has been observed in protostellar jets (e.g. Heathcote 
et al. 1996; Reipurth & Bally 2001). The knots in the mocked jet 
image move with the flow, consistent with observational findings 
(Eisloffel & Mundt 1992; Hartigan et al. 2005). This provides an 
alternative to the "internal shock" interpretation usually invoked 
to explain jet knots (e.g. Hartigan & Raymond 1993). 
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The magnetic structure of a jet undergoes a dramatic change 
if it becomes subject to violent kink instabilities. The formerly 
ordered helical structure is largely destroyed and the poloidal 
field becomes the dominating field component. However, mag- 
netic flux conservation still implies that the poloidal field de- 
clines faster with distance than the toroidal one. The toroidal 
field may thus, beyond the region covered by the here-presented 
simulations, become again dominant. This could, in principle, 
lead to a resurgence of instabilities until the Alfven speed has 
dropped below the sideways expansion speed of the jet (times 
a factor of order unity). From this point on the toroidal field is 
effectively frozen in the flow (cf. discussion in Paper 1). 

4.1. Collimation and jet environment 

Collimation of jets is popularly attributed to the "hoop stress" 
in the toroidal field component. This is misleading: though the 
stress contributed by the toroidal field can compress the config- 
uration near the axis (as observed in simulations), it eventually 
has to be taken up by an external agent (for a more extended dis- 
cussion, see Spruit 2008). In the simulations presented above, 
as well as in other work, this agent is an external medium sur- 
rounding the jet. It is included mostly because of limitations of 
the codes used, since the demands of conserving energy typi- 
cally cause numerical instabilities at low gas densities or low 
plasma-/?. 

The boundary between the jet and the external medium ac- 
tually expands due to pressure exerted by the toroidal field. The 
role of a (material) external medium in confining the jet can be 
taken over by a magnetic field, if it is able to transfer stress in the 
jet's toroidal field to the surface of the accretion disk. One might 
imagine that a toroidal field extending around the jet might serve 
this role. The high Alfven speeds in this field, however, would 
make it violently unstable to non-axisymmetric instabilities, as 
the early history of magnetic configurations for controlled fu- 
sion (linear pinches) testifies. This obstacle does not become 
evident in the axisymmetric models in the literature. A more 
realistic possibility, proposed by Shu et al. (1995), is that the 
disk's poloidal field, assumed for launching the jet in the inner 
regions, actually extends to much larger distances in the disk. 
Deformation of the field can take up the lateral stress exerted by 
the jet. The collimation of the jet would then be directly related 
to the properties of this poloidal field. High degrees of collima- 
tion would be most easily achieved in disks with a large ratio of 
outer to inner disk radius (Spruit et al. 1997). Numerical simu- 
lations at much lower plasma-/? than currently feasible would be 
needed to study this form of "poloidal collimation". 

As pointed out in the introduction, instabilities may take 
some distance to travel to become effective. The results pre- 
sented above show that a distance of the order of 1000 times 
the source size of the jet is needed to capture the dissipation of 
the toroidal field. I surmise this to be the reason why the effect of 
instabilities is not as noticeable in the works of Anderson et al. 
(2006) and Ouyed et al. (2003), rather than the lack of external 
confinement of these winds. 



"hot spot". These processes extract their energy directly from the 
bulk kinetic energy of the flow. The effect of internal instabili- 
ties deriving from the free energy in the toroidal field is much 
less destructive, since these mainly redistribute internal energy 
forms within the jet. The comparison between the 2.5D and the 
3D cases shows that the 3D jet widens by less than a factor 2 
as a result of internal instabilities (cf. Figs. 4 and 9). The non- 
axisymmetric nature of the instabilities thus does cause some 
interaction with the environment, but its consequences remain 
relatively benign. 

4.3. Cold flows 

Current 3D simulations are poorly equipped to handle flows in 
which the magnetic or kinetic energy density, or both, domi- 
nate over the thermal energy density. In magnetically dominated 
flows driven from actual accretion disks, the temperature of the 
plasma is often sufficiently low that magnetic energy density 
dominates over plasma pressure already at a short distance from 
the disk surface (Blandford & Payne 1982). As a consequence, 
the sonic point in the present results is further away from the 
source than would be the case in e.g. actual protostellar disks. 
It is possible that the onset of the instabilities in more strongly 
magnetically dominated flows would be faster, and their conse- 
quences even stronger than in the simulations presented here. 
Codes specially designed to handle such "cold" flows would be 
needed to verify this. 
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4.2. Disruption 

The possible presence of instabilities in jets sometimes raises 
concern about disruption. A complete dissipation of the jet into 
its surroundings is possible in principle through instabilities 
driven by the interaction of the jet with its surroundings, for ex- 
ample by Kelvin-Helmholtz instabilities, or its termination in a 



